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A  PRIMAL  DPG  METHOD  WITHOUT  A  FIRST  ORDER 

REFORMULATION 


L.  DEMKOWICZ  AND  J.  GOPALAKRISHNAN 
Communicated  by  Norbert  Heuer 

Abstract.  We  show  that  it  is  possible  to  apply  the  DPG  methodology  without  reformu¬ 
lating  a  second  order  boundary  value  problem  into  a  first  order  system,  by  considering  the 
simple  example  of  the  Poisson  equation.  The  result  is  a  new  weak  formulation  and  a  new 
DPG  method  for  the  Poisson  equation,  which  has  no  numerical  trace  variable,  but  has  a 
numerical  flux  approximation  on  the  element  interfaces,  in  addition  to  the  primal  interior 
variable. 


1.  Introduction 

A  typical  approach  to  derive  a  discontinuous  Petrov  Galerkin  (DPG)  method  for  a  second 
order  boundary  value  problem,  as  presented  in  [5,  6],  is  to  reformulate  the  problem  as  a  first 
order  system,  multiply  by  test  functions,  and  integrate  every  equation  by  parts.  In  this  note, 
we  show  that  it  is  not  necessary  to  reformulate  higher  order  equations  into  first  order  systems 
to  apply  the  DPG  methodology.  This  methodology  applied  to  a  simple  transport  equation 
was  presented  in  [4]  but  the  present  paper  uses  the  generalized  methodology  described  in  [5] . 

The  DPG  methodology  relies  on  automatically  computed  test  spaces.  To  describe  it, 
suppose  we  want  to  approximate  x  G  X  satisfying 

b{x,y)  =  l(y),  VyeY.  (l.l) 

Here  X  and  Y  are  real  Hilbert  spaces,  &(•,•)  :  X  x  Y  — >  M  is  a  continuous  bilinear  form, 
and  /(•)  :  Y  — >  M  is  a  continuous  linear  form.  We  use  a  closed  subspace  X h  C  X  for 
approximation  and  a  closed  subspace  Yr  C  Y  to  define  the  test  spaces  of  the  DPG  method. 
Typically  Xy  and  Yr  are  finite-dimensional  spaces  and  the  indices  h  and  r  are  related  to 
their  dimension.  Let  Tr  :  X  — y  Yr  be  defined  by  ( Trw,y)Y  =  b(w,y)  for  all  y  G  Yr.  Here 
and  throughout,  the  inner  product  (and  norm)  on  any  Hilbert  space  W  is  denoted  by  (-,  -)w 
(and  ||  •  || w,  respectively).  The  DPG  method  uses  the  test  space  Y£  =  Tr(Xh )  and  computes 
an  approximation  Xy  G  Xy  satisfying 

b(xh,y)  =  l(y)  VyeY,;.  (1.2) 

This  defines  a  Petrov-Galerkin  method  as  Xy  ^  Y£  in  general.  It  is  clear  from  this  abstract 
description  of  the  method  that  the  DPG  methodology  applies  very  generally  to  various 
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bilinear  forms.  In  particular,  there  is  no  need  to  rewrite  a  second  order  boundary  value 
problem  as  a  first  order  system  to  formulate  the  method  (1.2). 

Nonetheless,  for  the  method  to  be  practical,  the  operator  Tr  must  be  local.  This  is  usually 
ensured  by  finding  a  variational  setting  where  Yr  is  a  DG  space,  i.e.,  it  admits  functions  with 
no  continuity  constraints  across  mesh  element  interfaces.  For  the  Dirichlet  problem,  we  have 
shown  in  [6]  that  this  can  be  achieved  by  rewriting  it  as  a  first  order  system.  In  this  note,  we 
show  that  it  can  also  be  accomplished  without  using  the  first  order  reformulation.  In  doing 
so,  we  obtain  perhaps  the  simplest  example  illustrating  the  ideas  of  the  DPG  method. 

By  way  of  preliminaries,  we  recall  one  of  the  basic  theorems  used  to  analyze  DPG  methods, 
proved  in  [7].  It  uses  these  three  assumptions: 

(1)  Uniqueness: 

(to  6  I  :  b(w,v )  =  0  Vv  G  Y }  =  {0}.  (1.3a) 

(2)  Inf-sup  and  Continuity:  There  exist  C\ ,  C2  >  0  such  that 

CMy  <  sup  <  C2 |M|y.  (1.3b) 

wex  \\w\\x 

(3)  There  is  a  linear  operator  IJ  :  Y  — *  Yr  and  a  Cn  >  0  such  that  for  all  Wh  G  Xh  and 
all  v  G  Y, 

b(wh,v  -  ITv)  =  0,  ||ilu||y  <  Cn\\v\\Y-  (1.3c) 


Theorem  1.1  (see  [7]).  Suppose  Assumptions  (1.3a)  and  (1.3b)  hold.  Then  there  is  a  unique 
x  G  X  solving  (1.1),  i.e.,  Problem  (1.1)  is  well  posed.  If,  in  addition,  Assumption  (1.3c) 
holds,  then 


II  <r  C*Cn  •  f  II  II 
x-xh\\x<—^ —  lnf  \\x-w\\x- 

L-i  w€Xh 


(1.4) 


In  the  next  section,  we  present  the  new  weak  formulation  of  the  Dirichlet  problem.  In 
Section  3,  we  verify  the  first  two  assumptions  of  the  above  theorem,  thus  proving  that  the 
new  formulation  is  well-posed.  In  this  proof  we  use  that  a  piecewise  harmonic  function  can 
be  controlled  by  its  jumps  and  the  jumps  of  its  normal  derivatives,  a  fact  established  in  [6]. 
The  novelty  lies  in  the  bounds  of  gradients  of  piecewise  H 1  -functions  using  the  orthogonal 
projection  into  Vi?o(l2).  In  Section  4,  we  give  a  convergence  result  for  the  DPG  method 
applied  to  the  new  formulation.  We  conclude  in  Section  5  with  numerical  reports. 


2.  The  formulation 

Let  12  be  an  open  subset  of  M.N  (N  =  2  or  3)  with  Lipschitz  connected  boundary.  The 
boundary  value  problem  under  consideration  is 

—A u  =  f  on  i2,  (2.1a) 

u  —  0  on  812.  (2.1b) 

The  weak  formulation  of  this  problem  that  we  intend  to  study,  is  defined  using  a  mesh  12 h 
of  12.  We  assume  that  I2h  is  a  disjoint  partitioning  of  f2  into  open  elements  K  such  that  the 
union  of  their  closures  is  the  closure  of  12.  The  element  boundaries  d K  (for  all  the  elements 
K  G  12 h)  form  the  collection  dI2h-  We  assume  that  d K  is  Lipschitz  for  all  K  G  I2h,  so 
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that  we  may  use  trace  theorems  on  each  element,  but  the  shape  of  the  elements  is  otherwise 
arbitrary. 

The  new  Petrov-Galerkin  weak  formulation  of  this  problem  is  as  follows:  Find  (u,  qn)  G  X 
satisfying 

(V«,  Vu)i4  -  (qn,  v)anh  =  (/,  v)n  Vu  G  Y  (2.2) 

where  the  spaces  and  the  notations  used  are  defined  below. 

First,  define  a  broken  Sobolev  space  (where  as  usual  we  use  the  adjective  “broken”  for 
spaces  admitting  functions  with  interelement  discontinuities)  by  Ff1(i7/l)  =  {v  :  v\k  G 
if1  (if),  \/K  G  i?h},  and  normed  naturally  by 

M  m{Qh)  =  (v,  v)nh  +  (Vu,  Vv)nh.  (2.3) 

The  derivatives  above,  in  (2.2),  and  in  such  notations  throughout,  are  always  calculated 
element  by  element,  and 

(h«)ft  =  i^W)d(2h  =  w)l/2,dK- 

I<£fih  K&Qh 

where  denotes  the  L2(K)- inner  product  and  (£,-)i/2,di<  denotes  the  action  of  a  func¬ 

tional  i  in  if_1/2(i dK).  We  will  also  use  ||r||^  to  denote  the  norm  (r,  rfj^.  Next,  define  the 
space  of  numerical  fluxes  by 

H~l/2{dnh)  =  {ri  G  1^2((9 K)  :  3q  G  H( div,  Q)  such  that 

V\dK  =  Q  ■  n\ 9K,  V K  G  flh}i 

where  n  generically  denotes  the  unit  outward  normal  of  any  domain  under  consideration  -  for 
instance,  above  it  denotes  the  unit  outward  normal  on  the  boundary  of  a  mesh  element  K . 
This  space  is  normed  by 

||fn||^-i/2(5j7/i)  =  inf  {||g||//(div,D)  :  Vg  G  if  (div,  Q)  such  that  rn\dK  =  q  ■  n\dK}.  (2.4) 

With  these  notations,  we  can  now  complete  the  definition  of  the  weak  formulation  (2.2)  by 
setting  the  trial  and  test  spaces  in  (2.2)  by 

X  =  H(\(Q)  x  H-l/2(dQh) 

Y  =  H\nh). 

The  norm  in  Y  is  defined  by  (2.3)  and  the  norm  in  X  is  defined  by 

ll(w,rn)||x  =  l|Vw||  2L2{Q)  +  ||rn||^_i/2(aflh). 

This  completes  the  definition  of  all  the  notations  that  appeared  in  (2.2). 

At  this  point,  we  note  that  a  variational  problem  very  similar  to  (2.2)  can  be  found  in  [2, 
p.141] ,  by  the  name  of  a  primal  hybrid  method.  There,  it  is  used  as  motivation  to  introduce 
hybrid  and  non-conforming  methods.  An  important  difference  between  their  formulation 
and  (2.2)  is  that  while  theirs  is  standard  Galerkin  (or  Bubnov-Galerkin) ,  ours  is  a  Petrov- 
Galerkin  formulation  since  X  ^  Y .  We  now  proceed  to  prove  that  (2.2)  is  well-posed  in  the 
norms  introduced  above. 
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3.  Wellposedness 

Let  us  denote  the  bilinear  form  in  (2.2)  by 

b(  ( w,fn ),  v)  =  (Vw,  Vu)^  -  (■ rn,v)dnh 

for  any  ( w,  rn )  G  X  and  any  v  G  Y.  The  wellposedness  of  (2.2)  follows  from  Theorem  1.1, 
once  we  verify  (1.3a)  and  (1.3b). 

Lemma  3.1  (LIniqueness).  Any  ( w,rn )  G  X  satisfying 

b(  (w,rn),  v )  =  0  Vu  G  y  (3.1) 

vanishes. 


Proof.  Choosing  v  G  ‘D(iL)  on  any  hxed  mesh  element  K  G  Qh,  equation  (3.1)  implies  that 
— A(w|a-)  =  0.  Then  r  =  — Vw|a-  is  in  H( div,/l).  Therefore,  considering  a  general  v  eY, 
we  may  integrate  (3.1)  by  parts  element  by  element,  and  use  V  ■  (r \K)  =  0  to  obtain 

(r  •  n  -  r n,  =  0  Vu  G  iL1^). 


Choosing  v  supported  only  on  K,  this  shows  that  r  ■  n  —  rn  on  d K.  Furthermore,  for  any 
f  G  D(l7),  integrating  by  parts  element  by  element, 

-(c  V<f>)n  =  (V  •  r,  (f>)ah  +  (r  •  n,  f)dnh  =  (V  •  r,  +  (rn, 


Now,  by  definition  of  H  1/2(c)j7/l),  there  is  a  q  G  //(div,  17)  such  that  q  ■  n  —  rn  on  each 
d K.  Hence  integrating  by  parts,  first  element  by  element,  and  then  over  all  17,  we  have 
(' rn,4>)dnh  =  ( q-n,<t>)dnh  =  (V  ■  q,4>)n  +  (<lV0)«  =  {q-n,4>}en  =  0.  In  particular,  this 
proves  that  the  distributional  divergence  of  r  satisfies  (V  •  r){4>)  =  —  (r,  V0)j?  =  (Vt,  4>)ah, 
so  r  =  —Xiv  G  H(div,  17).  Therefore,  Aw  =  0  on  17  and  since  w  G  //(]  (17),  we  conclude  that 
vj  =  0.  Moreover  on  the  boundary  of  each  element,  rn  =  —dw/dn  =  0.  □ 


Next,  to  prove  (1.3b),  we  start  by  defining 

(r  •  n,  4>)anh 


[ t  ■  n 


def 

=  sup 


vn 


9Qh  ll^lli/qD) 

def  (f'm'V)oQh 


dQh 


(r-n,v)d{2h 

sup  - — r-n -  =  sup  rj  rj - . 

rn£H~1/2(dnh)  I  rn\\H-1/‘2(dQh)  r£H(div,0)  ||r  ||  H(div,D) 


=  sup 


(3.2a) 


(3.2b) 


The  last  equality  is  a  consequence  of  the  definition  of  the  quotient  norm  in  (2.4). 

An  important  ingredient  in  our  analysis  is  L2(l7)-orthogonal  projection  into  Vi7o(l7), 
which  we  denote  by  P: 

(• Pq ,  V0)fi  =  (q,  V0  G  Hl0{Q).  (3.3) 


Lemma  3.2.  There  is  a  positive  constant  C  independent  of  Qh  such  that  for  all  v  in  i71(l7/)  ), 

llull«k  <  C  (llPVyl|Dh  +  iHlsflJ  • 

Proof.  This  can  be  proved  along  the  same  lines  as  [6,  Lemma  4.2],  □ 


Lemma  3.3.  There  is  a  positive  constant  C  independent  of  Qh  such  that  for  all  v  in  H1(Qh), 

uv«lk  <  ll^ik+clHL,- 
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Proof.  By  the  definition  of  P,  there  is  z  G  H(\  (i?)  such  that  PWv  =  Vz.  Let  £  —  v  —  z  and 
define  r\K  =  —  We\k  on  all  K  G  17/j.  Then,  (3.3)  implies 

(r,V0)A-  =  O  (3.4) 

Choosing  <f>  G  D(/l),  we  immediately  find  that  V  •  (r| k)  —  0.  In  other  words, 

r  +  Vs  =  0  on  K , 

V  ■  r  =  0  on  K , 

Applying  a  known  result  [6,  Lemma  4.3]  for  locally  harmonic  functions,  we  thus  obtain 

IMk  =  l|Ve|k  <  (\ir  -n]\dnh  +  IM aah)  ■  (3-5) 

Now  we  claim  that  |  [r  •  n] \dQ  =  0.  This  is  because  we  may  integrate  by  parts  element  by 
element  to  conclude  from  (3.4)  that 

-(V  ■  r,  (f>)nh  +  (r  •  n,  4>)anh  =  (r  ■  n,  4>)dnh  =  0. 

Since  this  holds  for  all  (j)  G  H(\  (1?),  the  definition  (3.2a)  immediately  proves  the  claim.  It  also 
follows  from  (3.2b)  that  |  [s'lt]  | ;>Q  =  |  [vn]  \ dQ/  ■  With  these  facts  in  mind,  returning  to  (3.5), 
we  conclude  that 


IIVHk  <  +  ||V£||fil  <  ||PV«||fi,  +  c IMU, 

which  proves  the  lemma.  □ 

Lemma  3.4  (Inf-sup  condition).  There  exists  a  mesh-independent  C\  >  0  such  that 

b{  (w,rn),v) 


Ci||u||y  <  Slip 


Vu  G  Y. 


(3.6) 


Proof.  Let  s  denote  supremum  in  (3.6).  Choosing  rn  =  0  and  taking  the  supremum  over 
nontrivial  w  G  Po(i2)  we  obtain 

(Vw,  Vv)(2h 


S  >  sup 

weH'(n)  l|Vw||L2(fi) 


-  \\PVv\\v{{2). 


Next,  setting  w  =  0  and  maximizing  over  rn,  we  also  obtain 

-{rn,v)dnh 


s  > 


sup 


=  \vn 


dnh 


fneH-^idOh)  I \rn | \H- 1/2 (dnh) 
where  we  have  used  the  definition  (3.2b).  Thus,  we  have  shown  that 

2s  >  ||PVn||L2(r2)  +  |MLh 

for  all  v  G  Y.  The  required  inf-sup  condition  now  follows  from  Lemmas  3.2  and  3.3. 
Lemma  3.5  (Continuity).  There  is  a  mesh-independent  C2  >  0  such  that 
b(  (w,rn),v)  <  C2  ||(w,  rn) \\x  ||p||y,  V  (w,rn)  G  X,  v  G  Y. 


□ 


holds. 


(3.7) 
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Proof.  Since  it  is  immediate  from  the  definition  (3.2b)  that 

{?mv)df2h  5s  Il^n||_ff-1/2(gj7a)  |  \vn\ \gfih 

we  have 

b((w,rn),  v )  <  ||(if,rn)||x  (||Vu||^  +  ■ 

Hence,  it  suffices  to  prove  that 


|MLh<ClMI  y  VugF. 

But  using  the  definition  (3.2b)  again, 


TO 


(r  •  n,  v)dQh 


=  sup 


(r,  Vu)/jh  +  (V-r,u)flh 


h  rEH( div,i?)  ||^*||i/(div,/2)  r  EH (div,i?)  ||^||i/(div,i?) 

we  find  that  an  application  of  Cauchy-Schwarz  inequality  proves  (3.8). 


(3.8) 


□ 


4.  The  DPG  approximation 

In  this  section,  we  assume  that  Flh  consists  of  a  conforming  shape  regular  finite  element 
mesh  of  simplices  and  h  =  rnax^g^  (diarn  K).  We  turn  to  the  DPG  approximation  of 
(u,  qn)  G  X,  denoted  by  (uh,  qn,h )•  It  lies  in  Xh,  which  is  set  to 

Xh  =  {(w,b)  G  X  :  w\k  G  Pfc+i(A),  fn\gK  G  Pk(dK)},  (4.1) 

where  k  >  0  is  an  integer,  Pk(K )  denotes  the  space  of  polynomials  of  degree  at  most 
k  restricted  to  K  and  Pk(pK )  denotes  the  set  of  functions  on  d K  whose  values  on  each 
(N  —  1)  dimensional  subsimplcx  F  (face)  of  K  is  a  polynomial  of  degree  at  most  k  on  F. 
Next,  we  set 

Yr  =  {v  G  Y  :  v\K  G  Pr(A^)}. 

With  these  settings,  we  consider  the  DPG  method  (1.2)  and  its  solution  (uh,  qn,h),  and  prove 
the  following  result. 

Theorem  4.1.  Suppose  r  >  k  +  N.  Then,  there  is  a  constant  C  >  0  independent  of  h  such 
that 

\\u-uh\\Hi{n)  +  \\q- qnA\\H-i/2(dnh)  <C  inf  (||u  -  wh\\Hi{n)  +  \\q  -  fn>h\\H-i,2idn  )) . 

( Wh,rn,h)£Xh 

Proof.  We  apply  Theorem  1.1.  We  only  need  to  verify  (1.3c).  This  can  be  verified  element 
by  element.  Let  K  G  1?/, .  It  suffices  to  show  that  there  exists  a  bounded  linear  map 
PI  :  H\K )  -A  Pr(K)  such  that 

(V%,V(t)  -  nv))s 4  -  {fin, hi  V  -  nv)ds 4  =  0  (4.2) 

for  all  v  G  Hl(K),  wh  G  Pk{K)  and  qUih  G  Pk(dK). 

We  now  recall  the  result  of  [7,  Lemma  3.2],  which  asserts  that,  on  any  simplex  K,  there 
exists  a  C  independent  of  hx  =  diarn  (A")  and  a  bounded  linear  operator  np+x  '■  Hl(K)  — > 
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Pp+n(K)  such  that 

{np+Nv  -  v,  qp-i)K  =  0,  \/qp-i  G  Pp-i(K),  (4.3) 

(np+ NV  -  v,  nP)dK  =  o,  V/ip  G  Pp(dK),  (4.4) 

\\np+Nv\\Hi(K)  <  C\\v\\Hi{K),  \/v  G  Hl(K). 

Since  r  >  k  +  N,  we  may  set  77  =  IJk+N.  Then, 

(V%,V(d  -  IIv))nh  -  (qn,h,v  -  nv)aQh 

=  -  (A wh,  v  -  IIv)oh  ~  (qn,h  -  n  ■  Xwh,  v  -  nv)dQh  =  0, 
which  proves  (4.2),  and  concludes  the  proof.  □ 

Remark  4.1.  The  above  analysis  can  be  easily  generalized  to  a  mesh  of  hypercubes.  For 
example,  consider  the  case  of  squares  (TV  =  2)  and  let  QPtP  denote  the  (tensor  product) 
space  of  polynomials  of  degree  at  most  p  in  each  variable.  Then,  a  generalization  of  the 
construction  in  [7,  Lemma  3.2]  shows  that  on  any  square  K,  there  is  a  bounded  linear 
operator  77^2  :  Hl[K)  — >•  Qk+i,k+2{K)  such  that 

(nsk^2v  -  v,  qk)K  =  0,  Vqk  G  Qk,k( K), 

{nlq+2v  -  V,  Pk)dK  =  o,  V/Xfc  G  Pk(dK). 

(Details  of  this  construction  on  squares  are  explicitly  written  down  in  [3].)  Therefore,  pro¬ 
ceeding  as  in  the  proof  of  Theorem  4.1,  we  find  that  for  any  integer  k  >  1,  the  error  estimate 
of  the  theorem  continues  to  hold  on  square  meshes  with 

Xh  =  {{w,rn)  G  X  :  w\k  G  Qk,k( K),  rn\gK  G  Pk(dK)},  (4.5) 

=  {v  G  Y  :  v\k  e  Qk+2,k+2(K)}.  (4.6) 

Note  the  difference  in  the  degree  of  w  with  that  in  (4.1). 


Remark  4.2.  It  is  easy  to  see  that  the  condition  number  k  of  the  resulting  symmetric  positive 
definite  DPG  matrix  system  is  0(h~2)  on  quasiuniform  meshes.  It  is  proved  abstractly  in  [7, 
Remark  2.3]  that  any  DPG  method  has 


k  < 


Ai  ClCl 


n 


A0  Cl  ’ 

where  Ao  and  Ai  are  positive  constants  for  which  the  vector  of  coefficients  x  in  the  basis 
expansion  of  any  x  G  Xh  satisfies  A0||a;||^>  <  ||a;||^  <  Ai||T||^>-  Using  the  norm  equivalences 
in  the  proof  of  [7,  Theorem  3.7],  or  the  ones  in  [1],  we  can  easily  see  that  Ai/A0  <  Ch~ 2  for  our 
method.  The  remaining  constants  C'0,  C\ ,  Cjj  have  already  been  proved  to  be  independent 
of  h  for  the  formulation  under  consideration. 


5.  Numerical  illustration 

We  now  report  a  few  numerical  results.  We  consider  square  meshes  of  two  dimensional 
domains  and  set,  for  integers  k  >  1, 

Xh  =  {(w,rn)  G  X  :  w\ k  £  Qk,k(X),  rn\dK  £  Pk-i(dK)}, 

Yr  =  {v  G  Y  :  v\K  e  Qk+2,k+2(K)}. 
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(a)  The  square  case 
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L-shaped  domain:  h  convergence 
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(b)  The  case  of  the  L-shaped  domain 
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Figure  1.  Convergence  rates  on  a  square  and  L-shaped  domain 


The  above  defined  Xh  is  contained  in  (although  not  equal  to)  the  space  defined  in  (4.5),  so 
by  Remark  4.1,  the  error  estimate  of  Theorem  4.1  continues  to  hold.  We  choose  this  Xh  so 
that  the  best  approximation  rates  for  both  the  variables  are  comparable.  Convergence  rates 
in  terms  of  h  can  be  immediately  obtained  (the  best  approximation  error  can 

be  bounded  using  known  i/(div)-projections  of  the  flux  q  =  —Xu  as  in  [6,  Corollary  4.1]), 
namely, 

\\u  -uh]\m(n)  +  || q  -  Qn,h\\H-^(d{2h)  <  Chs  (\u\Hs+i{a)  +  |div  q\H°((2))  (5.1) 

for  all  s  <  k,  whenever  the  right  hand  side  norms  are  finite. 

We  illustrate  the  practical  manifestation  of  this  error  bound  by  solving  the  Dirichlet  prob¬ 
lem  (2.1)  on  the  unit  square  with  /  set  such  that  the  exact  solution  is  u  =  sin(7nc)  sin(7n/). 
Due  to  the  difficulties  in  computing  the  H~l'2(dQh)~ norm,  we  only  report  the  error  in  the 
iL1(i7)-norm  of  the  primal  variable: 

_  || U  —  Uh\\H1(fl) 

IMIiW«) 

The  computations  were  done  on  a  sequence  of  meshes  with  h  =  l/2n,  from  n  =  1, ...  ,6.  A 
plot  of  eh  versus  the  number  of  degrees  of  freedom  (which  is  0(/r~2))  is  shown  in  Figure  1(a). 
Clearly,  decreases  at  the  rates  predicted  by  (5.1)  for  each  k.  By  cross  connecting  the 
multiple  curves,  one  can  also  see  that  the  method  exhibits  exponential  decrease  of  error  with 
k  (i.e. ,  p-convergence),  even  though  that  is  not  covered  by  Theorem  4.1. 

Since  the  solution  in  above  experiment  was  infinitely  smooth,  the  maximal  rate  of  con¬ 
vergence  with  s  =  k  obtainable  in  (5.1)  was  observed  in  Figure  1(a).  In  contrast,  if  we 
approximate  the  singular  solution  u(r,6)  =  r2//3sin  (|(d  +  ^))  on  an  L-shaped  domain  us¬ 
ing  the  DPG  method,  then  as  seen  from  Figure  1(b),  due  to  the  regularity  limit,  we  only 
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observe  convergence  at  0(h 2/3).  Moreover,  the  convergence  rate  does  not  improve  even  if 
k  is  increased.  Since  the  solution  is  in  Hs+l((2)  for  s  <  2/3,  these  observations  are  also  in 
accordance  with  (5.1).  Note,  however,  that  by  cross  connecting  the  multiple  curves,  one  can 
observe  that  the  p-refinements  deliver  also  an  algebraic  convergence  with  a  doubled  rate  com¬ 
pared  with  /?, -refinements.  This  is  consistent  with  a  similar  behavior  of  standard  conforming 
elements. 
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